Автоматический подход к поиску неточностей

В любой работе с данными с ростом объемов информации, особенно в больших наборах данных, становится все сложнее находить выбросы и ошибки, которые могут исказить результаты анализа. Исследователи в первую очередь полагаются на визуальный осмотр данных, что удобно для небольших датасетов, однако плохо масштабируется с увеличением количества колонок и наблюдений.

С увеличением объема данных и их сложности стало очевидно, что нужны автоматизированные инструменты. В этом проекте мы рассматриваем наборы инструментов для автоматического поиска выбросов и ошибок в датасетах. Эти инструменты помогут ускорить анализ и сделать его более точным, позволяя исследователям сосредоточиться на более интересных задачах, сделав рутинное ковыряние в данных более эффективным и менее времязатратным.

В качестве примера данных загрузим датасет из проекта 1 (убедитесь, что загрузили эти данные в папку с этим Rmd файлом).

# main_dir <- dirname(rstudioapi::getSourceEditorContext()$path) 
# setwd(main_dir)
day_df <- read_csv("day.csv")
## Rows: 731 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (15): instant, season, yr, mnth, holiday, weekday, workingday, weathers...
## date  (1): dteday
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(day_df)

Рассматриваемые подходы и пакеты:

  1. dataMaid: Помогает автоматизировать проверку качества данных. Он может находить пропуски, ошибки в типах данных, выбросы, а также автоматически создавать отчёты по результатам анализа.
  2. skimr: Пакет для быстрого суммарного обзора данных, включая пропущенные значения, распределение переменных и типы данных.
  3. outliers: содержит тесты для оценки достоверности наличия выбросов.
  4. EnvStats: позволяет определить и оценить потенциальные выбросы.
  5. rstatix: позволяет искать многомерные выбросы на основе расстояния Махаланобиса и скорректированного расстояния Махаланобиса.
  • boxplot: встроенная функция, позволяющая строить графики и визуально оценивать наличие выбросов.

1. dataMaid

dataMaid простой пакет, который в 1 клик создаёт rmd-отчёт по вашему датафрейму и конвертирует его в pdf, docx, html.

Большим плюсом является кастомизируемость пакета, в него можно добавить свои функции для автоматического использования, что подробно описано в документации:

The customizability features of dataMaid essentially comes down to writing such summaryFunctions, visualFunctions and checkFunctions.

makeDataReport(day_df)
## The default of 'doScale' is FALSE now for stability;
##   set options(mc_doScale_quiet=TRUE) to suppress this (once per session) message
## Data report generation is finished. Please wait while your output file is being rendered.
# генерируем отчет в html
makeDataReport(day_df, output = "html", replace = TRUE)
## Data report generation is finished. Please wait while your output file is being rendered.
# подробно о функции
?makeDataReport()

Примеры работы:

Саммари по переменным
Саммари по переменным
Подробно о каждой переменной с возможными выбросами
Подробно о каждой переменной с возможными выбросами

Этот пакет удобен из-за нетрудоёмкости первичного анализа и наглядности генерируемых отчётов, однако не позволяет глубоко закопаться в предоставляемые данные.

2. Skimr

Пакет skimr предназначен для предоставления сводной статистики о переменных в дата фреймах, тибблах (tibble), таблицах данных (data tables) и векторах. В базовой библиотеке R наиболее похожую функцию выполняет summary() для векторов и дата фреймов.

Основной функцией пакета skimr является функция skim(), которая предназначена для работы с (группированными) датафреймами, она пытается преобразовать другие объекты в датафреймы, если это возможно. Подобно функции summary(), метод skim() для датафреймов отображает результаты для каждого столбца.

skim(day_df)
Data summary
Name day_df
Number of rows 731
Number of columns 16
_______________________
Column type frequency:
Date 1
numeric 15
________________________
Group variables None

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dteday 0 1 2011-01-01 2012-12-31 2012-01-01 731

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
instant 0 1.00 366.00 211.17 1.00 183.50 366.00 548.50 731.00 ▇▇▇▇▇
season 0 1.00 2.50 1.11 1.00 2.00 3.00 3.00 4.00 ▇▇▁▇▇
yr 0 1.00 0.50 0.50 0.00 0.00 1.00 1.00 1.00 ▇▁▁▁▇
mnth 0 1.00 6.52 3.45 1.00 4.00 7.00 10.00 12.00 ▇▅▅▅▇
holiday 0 1.00 0.03 0.17 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▁
weekday 0 1.00 3.00 2.00 0.00 1.00 3.00 5.00 6.00 ▇▃▃▃▇
workingday 0 1.00 0.68 0.47 0.00 0.00 1.00 1.00 1.00 ▃▁▁▁▇
weathersit 0 1.00 1.40 0.54 1.00 1.00 1.00 2.00 3.00 ▇▁▅▁▁
temp 2 1.00 0.50 0.18 0.06 0.34 0.50 0.66 0.86 ▂▇▇▇▅
atemp 0 1.00 0.47 0.16 0.08 0.34 0.49 0.61 0.84 ▂▇▆▇▂
hum 5 0.99 0.63 0.14 0.00 0.52 0.63 0.73 0.97 ▁▁▆▇▂
windspeed 1 1.00 0.19 0.08 0.02 0.13 0.18 0.23 0.51 ▃▇▅▁▁
casual 0 1.00 848.18 686.62 2.00 315.50 713.00 1096.00 3410.00 ▇▆▂▁▁
registered 1 1.00 3658.72 1559.81 20.00 2502.25 3664.50 4783.25 6946.00 ▂▅▇▅▃
cnt 0 1.00 4504.35 1937.21 22.00 3152.00 4548.00 5956.00 8714.00 ▂▅▇▅▃

Формат результатов, возвращаемых функцией skim(), представляет собой единый широкий датафрейм, который объединяет результаты анализа и включает в себя несколько атрибутов, а также два важных столбца метаданных:

skim_variable: содержит имена переменных из исходного набора данных. В нашем случае skim_variable — это колонка с именами переменных (Sepal.Length, Sepal.Width и т.д.).

skim_type: обозначает тип данных этих переменных (например, numeric, factor и т.д.). В нашем случае это Variable type: numeric и factor.

Эти столбцы являются ключевыми для объектов класса skim_df и играют важную роль в организации результатов. Если удалить эти столбцы, то объект потеряет свои уникальные свойства и будет преобразован в обычный объект типа tibble. Для проверки, пренадлежит ли объект к классу skim_df, используется функция is_skim_df().

В отличие от функции summary.data.frame(), которая хранит статистику в виде таблицы, объект skim_df имеет свои особенности. Это различие важно, поскольку skim_df можно использовать в пайпах, что облегчает дальнейшую работу с данными. Например, Вы можете выбрать все средние значения переменных или получить сводную статистику для конкретной переменной. Например, получим сводную статистику для столбца Petal.Length:

skim(day_df) %>%
  dplyr::filter(skim_variable == "hum") # поддерживается большинство функций из пакета dplyr
Data summary
Name day_df
Number of rows 731
Number of columns 16
_______________________
Column type frequency:
numeric 1
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
hum 5 0.99 0.63 0.14 0 0.52 0.63 0.73 0.97 ▁▁▆▇▂

Помимо датафреймов во второй версии пакета skimr функция skim() может анализировать векторы и матрицы путем преобразования их в датафрейм подобно тому как это делает функция as.data.frame():

# Обработка векторов
hum <- day_df$hum # выделяем вектор из данных.

skim(hum)
Data summary
Name hum
Number of rows 731
Number of columns 1
_______________________
Column type frequency:
numeric 1
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
data 5 0.99 0.63 0.14 0 0.52 0.63 0.73 0.97 ▁▁▆▇▂ 
all.equal(skim(hum), skim(as.data.frame(hum)))
## [1] "Attributes: < Component \"df_name\": 1 string mismatch >"
## [2] "Component \"skim_variable\": 1 string mismatch"
# Вывод означает, что объекты, которые мы сравнивали с помощью функции all.equal(), очень похожи, но отличаются в одном атрибуте, конкретно в компоненте df_name.
# Обработка матриц
day_matrix <- as.matrix(day_df)

skim(day_matrix) # аналогична функциям summary.matrix() и colMeans()
Data summary
Name day_matrix
Number of rows 731
Number of columns 16
_______________________
Column type frequency:
character 16
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
instant 0 1.00 3 3 0 731 0
dteday 0 1.00 10 10 0 731 0
season 0 1.00 1 1 0 4 0
yr 0 1.00 1 1 0 2 0
mnth 0 1.00 2 2 0 12 0
holiday 0 1.00 1 1 0 2 0
weekday 0 1.00 1 1 0 7 0
workingday 0 1.00 1 1 0 2 0
weathersit 0 1.00 1 1 0 3 0
temp 2 1.00 9 9 0 498 0
atemp 0 1.00 9 9 0 690 0
hum 5 0.99 8 8 0 591 0
windspeed 1 1.00 9 9 0 649 0
casual 0 1.00 4 4 0 606 0
registered 1 1.00 4 4 0 678 0
cnt 0 1.00 4 4 0 696 0

Пакет skimr так же предоставляет дополнительные функции для работы с данными.

Функция skim_tee() в пакете skimr позволяет получать ту же сводную информацию, что и skim(), но при этом возвращает исходный, неизменённый датафрейм. Это позволяет продолжить работу с оригинальными данными после получения сводной статистики. Посмотрим статистику для day_df, а затем отфильтруем датафрейм по рабочим дням:

day_working <- day_df %>%
  skim_tee() %>%
  dplyr::filter(workingday == "1")
## ── Data Summary ────────────────────────
##                            Values
## Name                       data  
## Number of rows             731   
## Number of columns          16    
## _______________________          
## Column type frequency:           
##   Date                     1     
##   numeric                  15    
## ________________________         
## Group variables            None  
## 
## ── Variable type: Date ─────────────────────────────────────────────────────────
##   skim_variable n_missing complete_rate min        max        median    
## 1 dteday                0             1 2011-01-01 2012-12-31 2012-01-01
##   n_unique
## 1      731
## 
## ── Variable type: numeric ──────────────────────────────────────────────────────
##    skim_variable n_missing complete_rate      mean        sd      p0      p25
##  1 instant               0         1      366       211.      1       184.   
##  2 season                0         1        2.50      1.11    1         2    
##  3 yr                    0         1        0.501     0.500   0         0    
##  4 mnth                  0         1        6.52      3.45    1         4    
##  5 holiday               0         1        0.0287    0.167   0         0    
##  6 weekday               0         1        3.00      2.00    0         1    
##  7 workingday            0         1        0.684     0.465   0         0    
##  8 weathersit            0         1        1.40      0.545   1         1    
##  9 temp                  2         0.997    0.496     0.183   0.0591    0.338
## 10 atemp                 0         1        0.474     0.163   0.0791    0.338
## 11 hum                   5         0.993    0.628     0.143   0         0.519
## 12 windspeed             1         0.999    0.191     0.0775  0.0224    0.135
## 13 casual                0         1      848.      687.      2       316.   
## 14 registered            1         0.999 3659.     1560.     20      2502.   
## 15 cnt                   0         1     4504.     1937.     22      3152    
##         p50      p75     p100 hist 
##  1  366      548.     731     ▇▇▇▇▇
##  2    3        3        4     ▇▇▁▇▇
##  3    1        1        1     ▇▁▁▁▇
##  4    7       10       12     ▇▅▅▅▇
##  5    0        0        1     ▇▁▁▁▁
##  6    3        5        6     ▇▃▃▃▇
##  7    1        1        1     ▃▁▁▁▇
##  8    1        2        3     ▇▁▅▁▁
##  9    0.5      0.656    0.862 ▂▇▇▇▅
## 10    0.487    0.609    0.841 ▂▇▆▇▂
## 11    0.628    0.730    0.972 ▁▁▆▇▂
## 12    0.181    0.233    0.507 ▃▇▅▁▁
## 13  713     1096     3410     ▇▆▂▁▁
## 14 3664.    4783.    6946     ▂▅▇▅▃
## 15 4548     5956     8714     ▂▅▇▅▃

Функция partition() разделяет результаты, полученные от skim(), на именованный список, где каждый элемент списка соответствует сводной информации для конкретного типа данных:

day_df %>%
  skim() %>%
  partition()

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dteday 0 1 2011-01-01 2012-12-31 2012-01-01 731

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
instant 0 1.00 366.00 211.17 1.00 183.50 366.00 548.50 731.00 ▇▇▇▇▇
season 0 1.00 2.50 1.11 1.00 2.00 3.00 3.00 4.00 ▇▇▁▇▇
yr 0 1.00 0.50 0.50 0.00 0.00 1.00 1.00 1.00 ▇▁▁▁▇
mnth 0 1.00 6.52 3.45 1.00 4.00 7.00 10.00 12.00 ▇▅▅▅▇
holiday 0 1.00 0.03 0.17 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▁
weekday 0 1.00 3.00 2.00 0.00 1.00 3.00 5.00 6.00 ▇▃▃▃▇
workingday 0 1.00 0.68 0.47 0.00 0.00 1.00 1.00 1.00 ▃▁▁▁▇
weathersit 0 1.00 1.40 0.54 1.00 1.00 1.00 2.00 3.00 ▇▁▅▁▁
temp 2 1.00 0.50 0.18 0.06 0.34 0.50 0.66 0.86 ▂▇▇▇▅
atemp 0 1.00 0.47 0.16 0.08 0.34 0.49 0.61 0.84 ▂▇▆▇▂
hum 5 0.99 0.63 0.14 0.00 0.52 0.63 0.73 0.97 ▁▁▆▇▂
windspeed 1 1.00 0.19 0.08 0.02 0.13 0.18 0.23 0.51 ▃▇▅▁▁
casual 0 1.00 848.18 686.62 2.00 315.50 713.00 1096.00 3410.00 ▇▆▂▁▁
registered 1 1.00 3658.72 1559.81 20.00 2502.25 3664.50 4783.25 6946.00 ▂▅▇▅▃
cnt 0 1.00 4504.35 1937.21 22.00 3152.00 4548.00 5956.00 8714.00 ▂▅▇▅▃

Функция yank() в свою очередь извлекает подтаблицу, соответствующую конкретному типу данных. Воспринимайте это как о dplyr::select для типов столбцов в исходных данных:

day_df %>%
  skim() %>%
  yank("numeric")

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
instant 0 1.00 366.00 211.17 1.00 183.50 366.00 548.50 731.00 ▇▇▇▇▇
season 0 1.00 2.50 1.11 1.00 2.00 3.00 3.00 4.00 ▇▇▁▇▇
yr 0 1.00 0.50 0.50 0.00 0.00 1.00 1.00 1.00 ▇▁▁▁▇
mnth 0 1.00 6.52 3.45 1.00 4.00 7.00 10.00 12.00 ▇▅▅▅▇
holiday 0 1.00 0.03 0.17 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▁
weekday 0 1.00 3.00 2.00 0.00 1.00 3.00 5.00 6.00 ▇▃▃▃▇
workingday 0 1.00 0.68 0.47 0.00 0.00 1.00 1.00 1.00 ▃▁▁▁▇
weathersit 0 1.00 1.40 0.54 1.00 1.00 1.00 2.00 3.00 ▇▁▅▁▁
temp 2 1.00 0.50 0.18 0.06 0.34 0.50 0.66 0.86 ▂▇▇▇▅
atemp 0 1.00 0.47 0.16 0.08 0.34 0.49 0.61 0.84 ▂▇▆▇▂
hum 5 0.99 0.63 0.14 0.00 0.52 0.63 0.73 0.97 ▁▁▆▇▂
windspeed 1 1.00 0.19 0.08 0.02 0.13 0.18 0.23 0.51 ▃▇▅▁▁
casual 0 1.00 848.18 686.62 2.00 315.50 713.00 1096.00 3410.00 ▇▆▂▁▁
registered 1 1.00 3658.72 1559.81 20.00 2502.25 3664.50 4783.25 6946.00 ▂▅▇▅▃
cnt 0 1.00 4504.35 1937.21 22.00 3152.00 4548.00 5956.00 8714.00 ▂▅▇▅▃

В документации представлены и другие функции из этого пакета, выше перечислены основные.

Отображение результатов skim()

Объект skim_df, получаемый от skim(), представлен в виде широкого датафрейма. По умолчанию его отображение осуществляется с помощью функции print.skim_df(). Для документов, создаваемых с использованием knitr, пакет skimr предоставляет специальный метод knit_print, который автоматически форматирует вывод. Чтобы использовать этот метод, последний оператор в вашем коде должен возвращать объект skim_df.

skim(day_df)
Data summary
Name day_df
Number of rows 731
Number of columns 16
_______________________
Column type frequency:
Date 1
numeric 15
________________________
Group variables None

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dteday 0 1 2011-01-01 2012-12-31 2012-01-01 731

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
instant 0 1.00 366.00 211.17 1.00 183.50 366.00 548.50 731.00 ▇▇▇▇▇
season 0 1.00 2.50 1.11 1.00 2.00 3.00 3.00 4.00 ▇▇▁▇▇
yr 0 1.00 0.50 0.50 0.00 0.00 1.00 1.00 1.00 ▇▁▁▁▇
mnth 0 1.00 6.52 3.45 1.00 4.00 7.00 10.00 12.00 ▇▅▅▅▇
holiday 0 1.00 0.03 0.17 0.00 0.00 0.00 0.00 1.00 ▇▁▁▁▁
weekday 0 1.00 3.00 2.00 0.00 1.00 3.00 5.00 6.00 ▇▃▃▃▇
workingday 0 1.00 0.68 0.47 0.00 0.00 1.00 1.00 1.00 ▃▁▁▁▇
weathersit 0 1.00 1.40 0.54 1.00 1.00 1.00 2.00 3.00 ▇▁▅▁▁
temp 2 1.00 0.50 0.18 0.06 0.34 0.50 0.66 0.86 ▂▇▇▇▅
atemp 0 1.00 0.47 0.16 0.08 0.34 0.49 0.61 0.84 ▂▇▆▇▂
hum 5 0.99 0.63 0.14 0.00 0.52 0.63 0.73 0.97 ▁▁▆▇▂
windspeed 1 1.00 0.19 0.08 0.02 0.13 0.18 0.23 0.51 ▃▇▅▁▁
casual 0 1.00 848.18 686.62 2.00 315.50 713.00 1096.00 3410.00 ▇▆▂▁▁
registered 1 1.00 3658.72 1559.81 20.00 2502.25 3664.50 4783.25 6946.00 ▂▅▇▅▃
cnt 0 1.00 4504.35 1937.21 22.00 3152.00 4548.00 5956.00 8714.00 ▂▅▇▅▃

3. Outliers

Пакет outliers содержит в себе несколько статистических тестов, позволяющих оценить наличие выбросов в принципе или задать их вручную и оценить достоверность этих выбросов.

  1. Тест Грабса

позволяет оценить, является ли максимальное или минимальное значение выбросом и оценивает уровень значимости. Работает на выборках размером больше \(6\).

По умолчанию считает нижний порог, чтобы посчитать статистику верхнего, нужно добавить параметр opposite = TRUE.

day_df = read.csv('day.csv')
grubbs.test(day_df$hum)
## 
##  Grubbs test for one outlier
## 
## data:  day_df$hum
## G = 4.39737, U = 0.97329, p-value = 0.003491
## alternative hypothesis: lowest value 0 is an outlier
grubbs.test(day_df$hum, opposite = TRUE)
## 
##  Grubbs test for one outlier
## 
## data:  day_df$hum
## G = 2.41300, U = 0.99196, p-value = 1
## alternative hypothesis: highest value 0.9725 is an outlier

В данном случае мы отвергаем, что значение \(0.9725\) это выброс, но принимаем, что \(0\) это выброс. Смотрим на графике:

out_day_df <- boxplot.stats(day_df$hum)$out
boxplot(day_df$hum)
mtext(paste("Outlier: ", paste(out_day_df, collapse = ", ")))

  1. Тест Диксона.

в отличие от теста Грабса, работает на маленьких выборках \(<=25\).

subday_df <- day_df[1:25, ] # отбираем из датасета первых 20 значений.
dixon.test(subday_df$hum)
## 
##  Dixon test for outliers
## 
## data:  subday_df$hum
## Q = 0.28209, p-value = 0.558
## alternative hypothesis: highest value 0.861667 is an outlier

В данном случае p.value слишком высок, значит мы отвергаем альтернативную гипотезу о том, что значение из нашей подвыборки является выбросом.

Можно перепроверить это с помощью boxplot.

out <- boxplot.stats(subday_df$hum)$out

boxplot(subday_df$hum)
mtext(paste("Outlier: ", paste(out, collapse = ", ")))

4. EnvStats

  1. тест Рознера

Пакет EnvStats содержит тест Рознера, который, в отличие от двух предыдущих тестов, используется для поиска сразу нескольких выбросов. Работает на выборках больше \(20\).

Здесь мы можем протестировать сразу несколько значений. Тест покажет, являются ли выбросы

Посмотрим на boxplot, чтобы определить потенциальное количество выбросов.

out_day_df <- boxplot.stats(day_df$hum)$out
boxplot(day_df$hum)
mtext(paste("Outlier: ", paste(out_day_df, collapse = ", ")))

При применении теста мы видим, что реальным выбросом является только 0. Второе значение не является выбросом. Мы можем узнать номер выброса в таблице при выводе теста, а вывод all.stats показывает нам произведённые расчёты. Функция автоматически удаляет NA из данных перед подсчётом квантилей.

test <- rosnerTest(day_df$hum, k = 2)
## Warning in rosnerTest(day_df$hum, k = 2): 5 observations with NA/NaN/Inf in 'x'
## removed.
test$all.stats

5. Rstatix

Пакет rstatix содержит функции для выявления одномерных и многомерных выбросов. Ищет с помощью boxplot и квартилей. Определяет выбросы и экстремальные выбросы с помощью квантилей (выше Q3 + 3 x IQR или ниже Q1 - 3 x IQR).

out_rstatix = identify_outliers(data = day_df, variable = 'hum')
out_rstatix

6. MVN

Пакет MVN ищет многомерные выбросы на основе расстояния Махаланобиса и скорректированного расстояния Махаланобиса.

С помощью функции mvn из пакета можно нарисовать графики Q-Q plot с выбросами на них по каждому из наблюдений, а также скорректированный график \(\chi\)-square.

# install.packages("MVN")
library(MVN)

sub_days = day_df[c(7, 10, 11, 12)]

result = mvn(data = sub_days, subset = "weekday", mvnTest = "hz",
             univariateTest = "AD", univariatePlot = "histogram",
             multivariatePlot = "qq", multivariateOutlierMethod = "adj",
             showOutliers = TRUE, showNewData = TRUE)

result$multivariateOutliers
## $`0`
##     Observation Mahalanobis Distance Outlier
## 408         408               95.516    TRUE
## 9             9               82.867    TRUE
## 23           23               49.306    TRUE
## 730         730               45.222    TRUE
## 380         380               41.246    TRUE
## 387         387               33.693    TRUE
## 429         429               27.186    TRUE
## 555         555               25.030    TRUE
## 415         415               24.226    TRUE
## 394         394               24.115    TRUE
## 205         205               23.205    TRUE
## 422         422               22.732    TRUE
## 51           51               18.688    TRUE
## 86           86               17.486    TRUE
## 16           16               16.551    TRUE
## 219         219               14.922    TRUE
## 268         268               14.211    TRUE
## 401         401               12.809    TRUE
## 352         352               10.556    TRUE
## 
## $`1`
##     Observation Mahalanobis Distance Outlier
## 10           10               82.805    TRUE
## 66           66               81.149    TRUE
## 367         367               75.787    TRUE
## 17           17               66.147    TRUE
## 24           24               62.073    TRUE
## 381         381               60.689    TRUE
## 87           87               54.167    TRUE
## 31           31               52.572    TRUE
## 52           52               50.578    TRUE
## 395         395               46.695    TRUE
## 430         430               44.501    TRUE
## 479         479               44.002    TRUE
## 409         409               42.973    TRUE
## 416         416               38.610    TRUE
## 731         731               35.834    TRUE
## 675         675               32.887    TRUE
## 206         206               32.539    TRUE
## 360         360               24.724    TRUE
## 192         192               21.450    TRUE
## 353         353               20.927    TRUE
## 45           45               18.502    TRUE
## 423         423               16.820    TRUE
## 388         388               14.583    TRUE
## 563         563               13.999    TRUE
## 374         374               13.338    TRUE
## 80           80               11.888    TRUE
## 269         269               11.480    TRUE
## 73           73               11.414    TRUE
## 402         402               10.818    TRUE
## 549         549                9.983    TRUE
## 
## $`2`
##     Observation Mahalanobis Distance Outlier
## 368         368               46.817    TRUE
## 39           39               28.605    TRUE
## 46           46               20.559    TRUE
## 88           88               15.249    TRUE
## 452         452               13.600    TRUE
## 53           53               12.317    TRUE
## 200         200               11.144    TRUE
## 
## $`3`
##     Observation Mahalanobis Distance Outlier
## 383         383               30.201    TRUE
## 369         369               29.074    TRUE
## 726         726               28.533    TRUE
## 26           26               24.137    TRUE
## 40           40               23.145    TRUE
## 677         677               21.085    TRUE
## 201         201               20.319    TRUE
## 362         362               20.022    TRUE
## 5             5               16.828    TRUE
## 208         208               14.637    TRUE
## 61           61               13.635    TRUE
## 222         222               12.921    TRUE
## 684         684               11.593    TRUE
## 698         698               11.515    TRUE
## 614         614               11.344    TRUE
## 33           33               11.211    TRUE
## 334         334                9.887    TRUE
## 54           54                9.813    TRUE
## 
## $`4`
##     Observation Mahalanobis Distance Outlier
## 202         202              110.358    TRUE
## 13           13               48.368    TRUE
## 727         727               47.870    TRUE
## 251         251               46.992    TRUE
## 265         265               37.518    TRUE
## 34           34               34.559    TRUE
## 69           69               29.710    TRUE
## 90           90               28.873    TRUE
## 41           41               25.555    TRUE
## 83           83               24.289    TRUE
## 62           62               22.232    TRUE
## 384         384               21.661    TRUE
## 342         342               17.089    TRUE
## 153         153               16.357    TRUE
## 20           20               15.495    TRUE
## 321         321               13.871    TRUE
## 405         405               11.745    TRUE
## 55           55               11.434    TRUE
## 545         545               11.323    TRUE
## 678         678               10.811    TRUE
## 
## $`5`
##     Observation Mahalanobis Distance Outlier
## 595         595             1914.826    TRUE
## 266         266               45.125    TRUE
## 203         203               40.278    TRUE
## 21           21               21.654    TRUE
## 252         252               19.709    TRUE
## 378         378               13.952    TRUE
## 
## $`6`
##     Observation Mahalanobis Distance Outlier
## 302         302               32.619    TRUE
## 421         421               32.106    TRUE
## 722         722               29.848    TRUE
## 694         694               27.903    TRUE
## 407         407               20.640    TRUE
## 8             8               19.133    TRUE
## 22           22               18.298    TRUE
## 386         386               17.634    TRUE
## 379         379               13.414    TRUE
## 204         204               12.757    TRUE
## 85           85               11.826    TRUE
## 435         435               11.407    TRUE
## 540         540               11.367    TRUE
## 351         351               11.349    TRUE
## 50           50                9.929    TRUE

* boxplot

Попробуем “полуавтоматический” поиск выбросов. Такой способ позволяет контролировать каждый этап. Воспользуемся встроенным методом boxplot.

С помощью Boxplot можно не только увидеть выбросы графически, но и сохранить их для последующего отсечения из данных. Для визуализации возьмём рандомные данные.

# создаём данные: с выбросами и без.
dirty_data <- c(4.358015, 4.489513, 4.560919, 4.613810, 4.599738, 4.621614, 4.633119, 4.616862, 4.754681, 4.849953, 4.945791, 5.019631, 4.805033, 4.989170, 5.024305, 5.065325, 4.970247, 4.998086, 5.096887, 4.977657, 4.888269, 3.479053, 2.878145)

data <- c(2.633213, 2.654674, 2.746650, 2.657763, 2.525229, 2.549804, 2.537088, 1.974909, 1.838017, 1.791683, 1.782088, 1.664908, 1.689402, 1.688826, 1.661763, 1.734322, 1.744875, 1.710471, 1.735690, 1.800677, 1.607354, 1.896810, 2.294757)

plot(data, dirty_data)

boxplot(data)

boxplot(dirty_data)

Визуально можем найти выбросы на графике. Чтобы удалить их, можно воспользоваться

boxplot не только рисует картинку, но и сохраняет все ее параметры в объекте, из которого мы можем их достать. Выбросы хранятся тут:

boxplot.stats(y)$out

Можно найти индексы точек выбросов.

ind <- which(dirty_data %in% boxplot.stats(dirty_data)$out)

ind
## [1] 22 23

Сохраним координаты точек выбросов в отдельном dataframe

out <- data.frame(data=data[ind], dirty_data=dirty_data[ind])
out

Если создать переменную с графиком (одним или несколькими), при выводе переменной мы получим следующие данные:

  • статистика: каждый столбец представляет нижний ус, первый квартиль, медиану, третий квартиль и верхний ус каждой группы.
  • n: количество наблюдений для каждой группы.
  • conf: каждый столбец представляет нижние и верхние крайние значения доверительного интервала медианы.
  • out: общее количество выбросов.
  • группа: общее количество групп.
  • имена: имена каждой группы.
res <- boxplot(dirty_data)

res
## $stats
##          [,1]
## [1,] 4.358015
## [2,] 4.606774
## [3,] 4.805033
## [4,] 4.983413
## [5,] 5.096887
## 
## $n
## [1] 23
## 
## $conf
##          [,1]
## [1,] 4.680948
## [2,] 4.929118
## 
## $out
## [1] 3.479053 2.878145
## 
## $group
## [1] 1 1
## 
## $names
## [1] "1"

И с помощью созданной переменной можно снова создать boxplot.

bxp(res)

Теперь визуализируем их на обычном графике:

plot(data, dirty_data, col='green', pch=18, ylim=c(0,max(dirty_data)))
points(out$data, out$dirty_data, col='red', pch=18)

Удалим точки из данных.

clean_data <- dirty_data[-ind]
data_data <- data[-ind]
clean_data
##  [1] 4.358015 4.489513 4.560919 4.613810 4.599738 4.621614 4.633119 4.616862
##  [9] 4.754681 4.849953 4.945791 5.019631 4.805033 4.989170 5.024305 5.065325
## [17] 4.970247 4.998086 5.096887 4.977657 4.888269
data_data
##  [1] 2.633213 2.654674 2.746650 2.657763 2.525229 2.549804 2.537088 1.974909
##  [9] 1.838017 1.791683 1.782088 1.664908 1.689402 1.688826 1.661763 1.734322
## [17] 1.744875 1.710471 1.735690 1.800677 1.607354

Сравним боксплоты с выбросами и без них:

boxplot(dirty_data, data_data)

boxplot(clean_data, data_data)

Сравним графики с выбросами и без них:

plot(data, dirty_data, col='green', pch=18, ylim=c(0,max(dirty_data)))
points(out$data, out$dirty_data, col='red', pch=18)

plot(data_data, clean_data, col='green', pch=18, ylim = c(0, max(clean_data)))